logo

Όλα τα απαραίτητα αρχεία βρίσκονται στο github.

1 Πριν ξεκινήσουμε

Βεβαιώστε ότι το όνομα του project σας, όπως και το file path ΔΕΝ περιέχουν:
1. Ελληνικούς χαρακτήρες
2. κενά (αντικαταστήστε τα κενά με _)

Επί παραδείγματι το file path σας θα μπορούσε να είναι: C:/Environmental_Data_Analysis/Geographical_Data.

3 Σχέση έκτασης - αριθμού ειδών

Η σχέση αυτή θεωρείται ένας από τους ελάχιστους νόμους της οικολογίας, καθώς έχει ισχύ σε όλες τις ταξινομικές βαθμίδες και όλους τους οργανισμούς (Triantis et al., 2012). Με πολύ απλά λόγια, θεωρεί ότι όσο αυξάνεται η έκταση μιας περιοχής, τόσο περισσότερα taxa αυτή φιλοξενεί.

Όμως υπάρχουν και άλλοι παράγοντες (κλιματολογικοί, τοπογραφικοί, γεωλογικοί και εδαφικοί), οι οποίοι μπορούν να επηρεάσουν τα πρότυπα ποικιλότητας και κατανομής των ειδών στον χώρο και τον χρόνο.

4 Περιοχή μελέτης

Σήμερα θα ασχοληθούμε με τα Ryukyus Islands, τα οποία είναι περισσότερο γνωστά ως Nansei Islands (στα Ιαπωνικά: Nansei-shoto, δηλαδή “Νοτιοδυτικά Νησιά”), ανήκουν στην Ιαπωνία. Αποτελούν ένα νησιωτικό τόξο, το οποίο βρίσκεται σε γεωγραφικό πλάτος 24ο - 31.0οΝ και εκτείνεται σε μήκος 1300 km σε έναν άξονα ΝΔ-ΒΑ, ξεκινώντας από το Kyushu και καταλήγει στην Taiwan. Το νοτιότερο νησί είναι το Yonaguni, ενώ το μεγαλύτερο η Okinawa. Τα μεγαλύτερα σε μέγεθος νησιά είναι κυρίως ορεινά, ενώ τα μικρότερα είναι κοραλλιογενή. Καθώς εντοπίζεται στην ζώνη μετάβασης από το τροπικό στο θερμό εύκρατο κλίμα (δηλαδή, βρίσκεται στην υποτροπική ζώνη), οι κλιματικές συνθήκες είναι ήπιες στο αρχιπέλαγος αυτό, καθ’όλη την διάρκεια του έτους, με μέση θερμοκρασία 15οC το χειμώνα και 28οC το καλοκαίρι. Βέβαια, υπάρχουν διαβαθμίσεις: στο Β τμήμα του αρχιπελάγους, το κλίμα χαρακτηρίζεται ως υγρό υποτροπικό, ενώ στο Ν τμήμα του ως τροπικό (ταξινόμηση κατά Koppen). Εξαιτίας της γεωγραφικής τους θέσης, στα νησιά αυτά δεν υπάρχει ξηρή περίοδος και τα περισσότερα εξ αυτών δέχονται περισσότερα από 2000 m βροχής, λόγω των μουσώνων και των καλοκαιρινών τυφώνων. Ως εκ τούτου, στα Ryukyus Islands δεν εμφανίζεται η υποτροπική ξηρή ζώνη, η οποία διαχωρίζει τα εύκρατα δάση από τις υγρές τροπικές περιοχές σε άλλα μέρη του πλανήτη και η φυσική τους βλάστηση αποτελείται από πλατύφυλλα αειθαλή δαση (Good, 1845; Maekawa, 1974; Kira, 1991).

Γεωλογικά, το αρχιπέλαγος διαιρείται σε 3 τμήματα:
1. Northern Ryukyus (τα νησιά βορείως του Tokara Strait)
2. Central Ryukyus (τα νησιά μεταξύ του Tokara Strait και του Kerama Gap)
3. Southern Ryukyus (τα νησιά νοτίως του Kerama Gap)

Εικόνα 1. Γεωτεκτονική διαίρεση του αρχιπελάγους Ryukyus.

Εικόνα 1. Γεωτεκτονική διαίρεση του αρχιπελάγους Ryukyus.

Τα Ryukyus είναι ηπειρωτικά νησιά και αποτελούν τις κορυφές της κορδιλλιέρας του αρχιπελάγους Ryukyus και έχουν δημιουργηθεί ως αποτέλεσμα της καταβύθισης της τεκτονικής πλάκας των Φιλιππίνων (Ota, 1998).

5 Εισαγωγή δεδομένων

Θα χρησιμοποιήσουμε τα ελεύθερα προσβάσιμα δεδομένα από τους Nakamura et al. (2009) για ορισμένα από τα νησιά αυτά. Τώρα ας εισαγάγουμε το αρχείο το οποίο περιέχει τα δεδομένα τα οποία μας ενδιαφέρουν να αναλύσουμε:

Ας δούμε πόσες στήλες και σειρές έχει:

## [1]   26 1815

Και τι είδους δεδομένα περιέχει:

## 'data.frame':    26 obs. of  1815 variables:
##  $ Codonacanthus_pauciflorus               : int  0 0 0 1 1 0 0 1 0 0 ...
##  $ Dicliptera_chinensis                    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Hemigraphis_reptans                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Hygrophila_salicifolia                  : int  0 0 0 1 0 0 0 1 0 0 ...
##  $ Justicia_procumbens                     : int  0 0 1 1 1 1 1 1 0 1 ...
##  $ Lepidagathis_formosensis                : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Lepidagathis_inaequalis                 : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Staurogyne_concinnula                   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Strobilanthes_glandulifer               : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Strobilanthes_tashiroi                  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Acer_capillipes                         : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Acer_insulare                           : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Acer_oblongum                           : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Acer_rufinerve                          : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Acer_sieboldianum                       : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Actinidia_hypoleuca                     : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Actinidia_polygama                      : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Actinidia_rufa                          : int  0 0 1 1 1 1 1 1 0 0 ...
##  $ Saurauria_roxburghii                    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Mollugo_pentaphylla                     : int  0 0 1 1 1 0 0 0 0 1 ...
##  $ Tetragonia_tetragonioides               : int  0 1 1 1 1 1 1 1 1 1 ...
##  $ Alangium_premnifolium                   : int  0 0 1 1 1 0 0 0 0 0 ...
##  $ Alisma_canaliculatum                    : int  0 0 0 1 0 0 0 1 0 1 ...
##  $ Caldesia_reniformis                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Sagittaria_aginashi                     : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Sagittaria_pygmaea                      : int  0 0 0 1 1 0 0 0 0 0 ...
##  $ Sagittaria_trifolia                     : int  0 0 0 1 1 1 0 0 0 0 ...
##  $ Achyranthes_aspera                      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Achyranthes_bidentata                   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Achyranthes_longifolia                  : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Deeringia_polysperma                    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Philoxerus_wrightii                     : int  0 0 0 1 1 0 0 1 0 0 ...
##  $ Crinum_asiaticum                        : int  1 0 1 1 1 1 1 1 1 1 ...
##  $ Lycoris_traubii                         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_crenata                            : int  0 0 1 1 1 0 1 1 0 0 ...
##  $ Rhus_ambigua                            : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Rhus_javanica                           : int  0 1 1 1 1 1 0 0 0 0 ...
##  $ Rhus_succedanea                         : int  0 0 1 1 1 1 0 1 0 0 ...
##  $ Polyalthia_liukiuensis                  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Anodendron_affine                       : int  0 1 1 1 1 1 0 0 0 0 ...
##  $ Cerbera_manghas                         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ecdysanthera_utilis                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Neisosperma_oppositifolia               : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Parsonsia_laevigata                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Trachelospermum_asiaticum               : int  0 1 1 1 1 1 1 1 1 1 ...
##  $ Trachelospermum_jasminoides             : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Ilex_dimorphophylla                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_ficoidea                           : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_goshiensis                         : int  0 0 0 1 1 0 0 0 0 0 ...
##  $ Ilex_hayataiana                         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_integra                            : int  1 1 1 1 1 1 1 1 1 0 ...
##  $ Ilex_kusanoi                            : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_liukiuensis                        : int  0 0 0 1 1 0 0 0 0 0 ...
##  $ Ilex_macrocarpa                         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_macropoda                          : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Ilex_maximowicziana                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ilex_pedunculosa                        : int  0 0 0 1 1 0 0 0 0 0 ...
##  $ Ilex_purpurea                           : int  0 0 0 1 0 0 0 0 0 0 ...
##  $ Ilex_rotunda                            : int  0 0 1 1 1 1 1 1 0 0 ...
##  $ Acorus_gramineus                        : int  0 0 1 1 1 0 0 0 0 0 ...
##  $ Alocasia_cucullata                      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Alocasia_gageana                        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Alocasia_odora                          : int  0 0 1 1 1 1 1 1 0 1 ...
##  $ Amorphophallus_kiusianus                : int  0 0 0 1 1 0 0 1 0 0 ...
##  $ Arisaema_heterocephalum                 : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Arisaema_japonicum                      : int  0 0 1 1 1 0 0 1 0 0 ...
##  $ Arisaema_kawashimae                     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Arisaema_longipedunculatum              : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Arisaema_ringens                        : int  0 0 1 1 1 1 1 1 1 1 ...
##  $ Arisaema_robustum                       : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Arisaema_sazensoo                       : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Arisaema_thunbergii                     : int  0 0 1 1 1 1 1 1 0 0 ...
##  $ Arisaema_yakualpinum                    : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Pinellia_ternata                        : int  0 0 1 1 1 0 0 0 0 0 ...
##  $ Pinellia_tripartia                      : int  0 0 0 0 1 0 0 0 1 1 ...
##  $ Rhaphidophora_liukiuensis               : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Rhaphidophora_pinnata                   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Typhonium_divaricatum                   : int  0 0 0 1 0 0 0 0 0 0 ...
##  $ Aralia_elata                            : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Dendropanax_trifidus                    : int  0 1 1 1 1 1 1 1 1 0 ...
##  $ Fatsia_japonica                         : int  1 1 1 1 1 1 1 1 1 0 ...
##  $ Hedera_rhombea                          : int  0 0 1 1 1 0 1 1 0 0 ...
##  $ Kalopanax_pictus                        : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Kalopanax_septemlobus                   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Schefflera_octophylla                   : int  0 0 1 1 1 1 1 1 0 0 ...
##  $ Aristolochia_debilis                    : int  0 0 0 1 0 0 0 0 0 0 ...
##  $ Aristolochia_kaempferi                  : int  0 0 0 1 1 0 0 0 0 0 ...
##  $ Aristolochia_liukiuensis                : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Aristolochia_zollingeriana              : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_caudigerum                       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_celsum                           : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_dissitum                         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_fudsinoi                         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_gelasinum                        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_gusuk                            : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_hatsushimae                      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_hirsuitsepalum                   : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ Asarum_lutchuense                       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Asarum_monodoriflorum                   : int  0 0 0 0 0 0 0 0 0 0 ...
##   [list output truncated]

6 Μετατροπή μήτρας παρουσίας/απουσίας

Όπως θα διαπιστώσατε, τα δεδομένα αυτά αποτελούν μια μήτρα παρουσίας/απουσίας ειδών (presence/absence matrix). Τέτοιου είδους δεδομένα είναι πολύ χρήσιμα για αναλύσεις σχετικά με την β-ποικιλότητα, όμως σήμερα εμείς χρειαζόμαστε μόνο τον αριθμό των ειδών σε ένα συγκεκριμένο νησί και όχι ποια είδη βρίσκονται σε κάθε νησί.

Ένας τρόπος προκειμένου να μετατρέψουμε μια μήτρα παρουσίας/απουσία ειδών στον συνολικό αριθμό ειδών τα οποία απαντώνται σε κάθε νησί, είναι ο ακόλουθος:

Island SR
Takeshima 163
Ioujima 158
Kuroshima 426
Tanegashima 839
Yakushima 1025
Kuchinoerabujima 313
Kuchinoshima 336
Nakanoshima 494
Gajajima 234
Tairajima 241
Suwanosejima 267
Akusekijima 318
Kodakarajima 153
Takarajima 413
Yokoatejima 67
Amamioshima 863
Kikaijima 446
Tokunoshima 675
Okinoerabujima 568
Yoronjima 358
Okinawajima 1011
Kumejima 597
Miyakojima 552
Ishigakijima 877
Iriomotejima 887
Yonagunijima 590

7 Εισαγωγή γεωχωρικών δεδομένων

Καθώς έχουμε πλέον στη μορφή που θέλουμε τα βιοτικά μας δεδομένα, επόμενο μας βήμα είναι να συγκεντρώσουμε αβιοτικά δεδομένα για τα νησιά αυτά.

7.1 Διανυσματικά δεδομένα

Για τον σκοπό αυτό, θα χρειαστεί πρώτα να εισαγάγουμε τα διανυσματικά αρχεία (shapefiles) για τα νησιά που μας ενδιαφέρουν.

Θα χρειαστεί να αλλάξετε το file path σε Shapefiles.

7.2 Ψηφιδοπλέγματα

Ας είσαγουμε τώρα και τα ψηφιδοπλέγματα τα οποία περιέχουν τις αβιοτικές μας μεταβλητές.

Αναλόγως των δυνατοτήτων του Η/Υ σας, υπάρχει περίπτωση οι ακόλουθες εντολές να διαρκέσουν μεγάλο χρονικό διάστημα. Μπορείτε να εισαγάγετε το αποτέλεσμα αυτών με την εντολή preds <- readRDS('Data/RDS/Ryukyus abiotic data.rds') & altitude <- readRDS('Data/RDS/Ryukyus altitudinal data.rds')

Ας μεταφορτώσουμε και υψομετρικά δεδομένα για τα νησιά αυτά:

8 Εξαγωγή γεωχωρικών δεδομένων

Τώρα μπορούμε να εξάγουμε τα αβιοτικά αυτά δεδομένα για κάθε νησί που εμπεριέχεται στην ανάλυση μας:

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |======================================================================| 100%
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |======================================================================| 100%

9 Αντιμετώπιση προβλημάτων

Εάν παρατηρήσετε τα δεδομένα για τα νησιά που μας ενδιαφέρουν, θα δείτε2, ότι για ορισμένα νησιά δεν έχουμε καθόλου (NA) δεδομένα. Αυτό μπορεί να οφείλεται είτε στο γεγονός ότι τα νησιά αυτά μπορεί να είναι πολύ μικρά, οπότε να μην υπάρχουν διαθέσιμα δεδομένα στο μέγεθος των ψηφιδοπλεγμάτων τα οποία χρησιμοποιούμε (resolution - ακρίβεια), είτε εξαιτίας του γεγονότος ότι χρησιμοποιούμε σχετικά αδρά ψηφιδοπλέγματα να μην έχουμε δεδομένα για αυτά τα νησιά, ανεξαρτήτως του μεγέθους τους3.

Υπάρχουν δύο αποδεκτές λύσεις:

  • Μετασχηματισμός των δεδομένων μας
  • Αφαίρεση των προβληματικών δεδομένων

Και οι δύο προσεγγίσεις έχουν τα πλεονεκτήματα και τα μειονεκτήματα τους. Η πρώτη προσέγγιση είναι προτιμότερη όταν έχουμε πλήθος διαθέσιμων δεδομένων (>5000 δείγματα) και το ποσοστό των προβληματικών δεδομένων είναι μικρό (συνήθως < 1-2%). Η δεύτερη προσέγγιση θεωρείται γενικά ασφαλέστερη, διότι δεν δημιουργούμε εκ του μηδενός δεδομένα. Σήμερα θα ακολουθήσουμε την δεύτερη προσέγγιση, όμως θα δούμε πως μπορούμε να μετασχηματίσουμε τα δεδομένα μας.

9.2 Αφαίρεση δεδομένων

Ας δούμε τώρα πως μπορούμε να αφαιρέσουμε τα όποια προβληματικά δεδομένα μας:

##  [1] "Akusekijima"  "Gajajima"     "Ioujima"      "Kodakarajima" "Kuchinoshima"
##  [6] "Tairajima"    "Takarajima"   "Takeshima"    "Yokoatejima"  "Yonagunijima"
## ============================================================================##
## Create a symbol to filter out all the Islands having NA values
## ============================================================================##
"%!in%" <- function(x, y) !(x %in% y)
## ============================================================================##


## ============================================================================##
## Now create a temporary presence/absence matrix removing any problematic islands
## ============================================================================##
plants_new <- plants %>% mutate(Island = rownames(.)) %>% as_tibble() %>% dplyr::select(Island, 
    everything()) %>% pivot_longer(!Island, names_to = "Taxon", values_to = "PA") %>% 
    filter(Island %!in% setdiff(imputed_vars$Island, drop_vars$Island)) %>% pivot_wider(names_from = Taxon, 
    values_from = PA)
## ============================================================================##


## ============================================================================##
## See which taxa are left after removing the problematic islands
## ============================================================================##
keepers <- colSums(plants_new[, -1]) %>% enframe() %>% filter(value > 0)
## ============================================================================##


## ============================================================================##
## Finally create the new PA matrix
## ============================================================================##
plants_new_pa <- plants_new %>% pivot_longer(!Island, names_to = "Taxon", values_to = "PA") %>% 
    filter(Taxon %in% keepers$name) %>% pivot_wider(names_from = Taxon, values_from = PA)
## ============================================================================##

10 Κατασκευή φυλογενετικού δένδρου

Προκειμένου να μπορέσουμε να υπολογίσουμε την φυλογενετική α-ποικιλότητα για το αρχιπέλαγος της Ιαπωνίας, θα χρειαστεί να κατασκευάσουμε πρώτα το φυλογενετικό δένδρο των ειδών που απαρτίζουν τις νησιωτικές φυτοκοινότητες. Για τον σκοπό αυτό, θα χρησιμοποιήσουμε εντολές από τις βιβλιοθήκες tidyr, dplyr & rtrees. Η τελευταία βιβλιοθήκη περιλαμβάνει τα φυλογενετικά δένδρα για αρκετούς φυτικούς οργανισμούς, καθώς και για τους ιχθύες, τα πτηνά και τα θηλαστικά. Θα χρησιμοποιήσουμε μόνο τα taxa τα οποία απαντώνται στα 16 νησιά του αρχιπελάγους της Ιαπωνίας που δεν είχαν κάποιο πρόβλημα με τις αβιοτικές μεταβλητές. Τα taxa αυτά βρίσκονται στον πίνακα που ονομάσαμε keepers. Θα χρειαστεί να προσθέσουμε στον πίνακα αυτόν το όνομα του γένους για κάθε φυτικό taxon το οποίο απαντάται στα 16 αυτά νησιά των Ryukyus και στη συνέχεια, θα χρησιμοποιήσουμε αυτόν τον πίνακα για να κατασκευάσουμε το φυλογενετικό μας δένδρο. Ακόμα και εάν το φυλογενετικό μας δένδρο περιέχει πολυτομίες, είναι ασφαλές να το χρησιμοποιήσουμε για μακροοικολογικές αναλύσεις (Li et al., 2019).

Αναλόγως των δυνατοτήτων του Η/Υ σας, υπάρχει περίπτωση η ακόλουθη εντολή να διαρκέσει μεγάλο χρονικό διάστημα. Μπορείτε να εισαγάγετε το αποτέλεσμα αυτής με την εντολή: jpn_tree <- readRDS('Data/RDS/16 Ryukyus tree.rds')

Στη συνέχεια, μπορούμε να αναπαραστήσουμε γραφικά το φυλογενετικό δένδρο, δίνοντας διαφορετικό χρώμα σε κάθε γένος:

10.1 Υπολογισμός φυλογενετικής α-ποικιλότητας

Στο σημείο αυτό, θα χρησιμοποιήσουμε μια εντολή από την βιβλιοθήκη picante, προκειμένου να κρατήσουμε μόνο εκείνα τα είδη τα οποία είναι κοινά μεταξύ του φυλογενετικού δένδρου και της μήτρας παρουσίας/απουσίας:

## [1] "Dropping taxa from the community because they are not present in the phylogeny:"
##  [1] "Saurauria_roxburghii"       "Carmona_retusa"            
##  [3] "Saionia_shinzatoi"          "Crataeva_falcata"          
##  [5] "Pseudostellaris_heterantha" "Dispanthus_uniflorus"      
##  [7] "Siegesbeckia_glabrescens"   "Siegesbeckia_orientalis"   
##  [9] "Cardamina_flexuosa"         "Pogonantherum_crinitum"    
## [11] "Spodipogon_cotulifer"       "Eusteralis_stellata"       
## [13] "Perillura_reptans"          "Akebi_quinata"             
## [15] "Dumbaria_villosa"           "Alectrorurus_yedoensis"    
## [17] "Caphalanthera_falcata"      "Heterozeuxine_nervosa"     
## [19] "Heterozeuxine_odorata"      "Luisaerides_liukiuensis"   
## [21] "Vexillabium_yakushimense"   "Hydyotis_verticillata"     
## [23] "phtheirospermum_japonicum"  "Cnidum_japonicum"          
## [1] "Dropping tips from the tree because they are not present in the community data:"
## [1] "Phtheirospermum_japonicum"

Τώρα μπορούμε να υπολογίσουμε την φυλογενετική α-ποικιλότητα, μέσω της βιβλιοθήκης PhyloMeasures:

11 Τελικό dataset

Τώρα πρέπει να κρατήσουμε μόνο εκείνα τα νησιά που δεν έχουν προβληματικά δεδομένα, να τα ενώσουμε με τα αβιοτικά και τα φυλογενετικά δεδομένα, καθώς και να υπολογίσουμε την έκταση του εκάστοτε νησιού:

##============================================================================##
## Create the final dataset
##============================================================================##
ryukyus_new <- ryukyus %>% 
  filter(Island %in% plants_new_pa$Island) %>% 
  full_join(., drop_vars) %>% 
  full_join(., ryukyus_shps) %>% 
  filter(Island %!in% setdiff(imputed_vars$Island, 
                              drop_vars$Island)) %>% 
  st_as_sf() %>% 
  
  ## Lambert Azimuthal Equal Area projection
  st_transform('+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs') %>% 
  
  ## Create the Area and PD variables
  mutate(Area = as.numeric(units::set_units(st_area(.), km^2)),
         PD = japan_pd,
         PD_P = japan_pd_p) %>% 
  
  ## Transform back to WGS84
  st_transform(4326) %>% 
  
  dplyr::select(Island, SR, PD, PD_P, Area, Elevation, everything()) %>% 
  
  ## It has zero-variance
  dplyr::select(-aridityIndexThornthwaite) %>% 
  
  ## Drop the geometry
  st_drop_geometry() %>% 
  
  ## Join these data with the Maximum and Minimum distance between each island
  ## pair
  inner_join(., 
             
             ## Calculate the Maximum distance between each island pair
             apply((st_distance(ryukyus_shps %>% 
                                     st_centroid(), 
                                   ryukyus_shps %>% 
                                     st_centroid())/1e+3) %>% 
                        as_tibble() %>% 
                        set_names(ryukyus_shps$Island) %>% 
                        units::drop_units() %>% 
                        na_if(., 0), 1, max, na.rm = T) %>% 
               enframe() %>% 
               rename(MaxD = value) %>% 
               mutate(Island = ryukyus_shps$Island) %>% 
               dplyr::select(Island, MaxD) %>% 
               
               full_join(., 
                         
                         ## Calculate the Minimum distance between each island
                         ## pair
                         apply((st_distance(ryukyus_shps %>% 
                                                 st_centroid(), 
                                               ryukyus_shps %>% 
                                                 st_centroid())/1e+3) %>% 
                                    as_tibble() %>% 
                                    set_names(ryukyus_shps$Island) %>% 
                                    units::drop_units() %>% 
                                    na_if(., 0), 1, min, na.rm = T) %>% 
                           enframe() %>% 
                           rename(MinD = value) %>% 
                           mutate(Island = ryukyus_shps$Island) %>% 
                           dplyr::select(Island, MinD)))
##============================================================================##

ryukyus_new
Island SR PD PD_P Area Elevation bio1 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19 bio2 bio3 bio4 bio5 bio6 bio7 bio8 bio9 annualPET climaticMoistureIndex continentality embergerQ growingDegDays0 growingDegDays5 maxTempColdest minTempWarmest monthCountByTemp10 PETColdestQuarter PETDriestQuarter PETseasonality PETWarmestQuarter PETWettestQuarter thermicityIndex MaxD MinD
Kuroshima 426 -1.8402810 0.9660340 16.30184 213.81479 176.0000 253.0000 102.00000 2636.000 516.0000 117.0000 48.00000 1079.0000 403.0000 717.0000 408.0000 68.00000 29.00000 5816.000 293.0000 59.00000 234.0000 221.0000 157.00000 1000.6500 0.6200000 16.90000 387.3800 76194.00 76194.00 125.0000 229.0000 10.000000 44.86000 58.52000 3070.150 113.6200 115.83000 360.0000 983.8313 34.73631
Tanegashima 839 -2.0300010 0.9830170 449.44967 89.10355 189.1470 264.5754 115.14228 2759.873 492.1473 108.4328 44.06544 1075.9171 374.1039 805.9562 374.1039 61.59619 26.50000 5791.570 302.0701 73.81043 228.0324 234.1318 115.14228 990.7673 0.6451378 16.82070 417.7372 81746.30 81746.30 135.3028 242.8109 11.500000 44.90803 44.90803 3037.074 113.2602 114.93804 398.3696 1040.1992 51.08020
Yakushima 1025 -0.2158943 0.5644356 508.57513 606.53490 158.5849 231.3554 86.39672 3318.087 572.3122 137.1453 41.59947 1270.1959 486.8870 964.3652 486.8870 61.16105 27.14583 5602.943 267.2355 46.68151 220.8119 199.3229 86.39672 909.2751 0.7249949 16.21991 520.0648 68752.66 68752.66 106.3139 210.1564 9.506935 40.78871 40.78871 2838.298 103.5724 105.73511 311.2719 989.5511 32.70489
Kuchinoerabujima 313 -1.6912781 0.9610390 37.17430 158.07548 184.0000 256.0000 114.00000 2812.000 543.0000 123.0000 45.00000 1105.0000 459.0000 724.0000 459.0000 64.00000 29.00000 5478.000 293.0000 74.00000 219.0000 196.0000 114.00000 1000.3900 0.6400000 15.80000 440.4100 79416.00 79416.00 135.0000 232.0000 12.000000 45.54000 45.54000 3059.810 113.1800 107.40000 393.0000 975.0262 32.70489
Nakanoshima 494 -1.0391167 0.8381618 35.72809 244.47767 185.0000 254.0000 118.00000 3042.000 564.0000 146.0000 43.00000 1199.0000 498.0000 748.0000 498.0000 61.00000 28.00000 5234.000 291.0000 80.00000 211.0000 194.0000 118.00000 978.1600 0.6800000 15.30000 494.1600 80010.00 80010.00 138.0000 233.0000 12.000000 45.69000 45.69000 2878.560 110.3400 103.21000 403.0000 906.5470 14.56124
Suwanosejima 267 -0.9998780 0.8271728 29.85171 191.07586 191.0000 260.0000 125.00000 3033.000 547.5000 148.5000 42.00000 1182.0000 499.5000 743.5000 499.5000 60.00000 28.00000 5202.000 295.5000 86.50000 209.0000 199.0000 125.00000 988.7250 0.6750000 15.10000 496.4800 82503.00 82503.00 144.5000 237.5000 12.000000 46.77500 46.77500 2862.030 111.2700 103.63500 422.0000 880.1082 18.14204
Amamioshima 863 -0.8816950 0.8081918 724.01156 168.94980 207.2345 269.9245 144.86514 2927.709 412.3831 145.5170 30.46525 964.3180 490.3122 837.2828 490.3122 57.82408 29.00000 4885.349 305.1555 109.35008 196.0141 219.1768 144.86514 1019.7703 0.6466562 13.94785 505.2497 89714.96 89714.96 165.8180 248.5766 12.000000 50.93358 50.93358 2752.129 112.8781 104.93128 482.6365 771.3206 53.54787
Kikaijima 446 -2.4427093 0.9930070 59.58188 34.49266 213.5000 276.5000 153.01189 2868.929 416.6492 143.9881 31.50000 959.9048 479.0000 802.7758 479.0000 58.50000 29.00000 4851.000 312.0000 115.50000 196.5000 226.0119 153.01189 1050.6402 0.6350000 13.84612 494.0876 93132.43 93132.43 174.0119 253.5000 12.000000 52.36262 52.36262 2846.930 116.6501 108.38107 506.0119 817.8176 53.54787
Tokunoshima 675 -1.0657777 0.8661339 253.48894 123.09407 213.5393 273.4286 153.53926 2409.559 319.2652 123.8041 27.50000 811.3153 411.0024 660.5954 411.0024 53.50000 28.00000 4683.245 305.9859 119.45537 186.5074 224.4286 153.53926 999.6211 0.5852746 13.27762 438.1796 92427.59 92427.59 173.5393 253.0934 12.000000 51.34734 51.34734 2598.477 110.1598 101.82562 506.3661 700.2143 55.69807
Okinoerabujima 568 -1.8609689 0.9640360 95.13222 72.94901 217.8094 276.4325 159.32750 2148.366 273.9207 111.8200 26.50000 729.3018 368.1002 580.7622 370.2275 51.50000 28.00000 4602.111 307.4325 125.80945 181.5000 228.7475 159.32750 994.5787 0.5364186 12.92500 401.9761 94279.97 94279.97 179.1342 255.8094 12.000000 52.08294 52.08294 2518.234 109.2814 100.45288 522.4284 648.5441 40.89212
Yoronjima 358 -2.2220711 0.9880120 20.92543 26.13938 222.0000 279.0000 163.00000 2139.000 274.0000 111.0000 28.00000 713.0000 362.0000 603.0000 362.0000 53.00000 28.00000 4556.000 311.0000 128.00000 183.0000 275.0000 163.00000 1026.5900 0.5200000 12.85000 396.0200 95922.00 95922.00 184.0000 258.0000 12.000000 54.10000 54.10000 2567.800 112.4200 115.28000 534.0000 616.5526 40.89212
Okinawajima 1011 -1.5138595 0.9350649 1236.14528 74.98082 220.2019 276.5807 161.14858 2216.502 290.1055 115.5279 31.38394 777.4536 367.2347 665.8982 367.2347 54.85507 29.50000 4513.869 309.0611 124.09077 183.2771 273.3103 161.14858 1038.3087 0.5243882 12.79125 404.2533 95252.11 95252.11 182.3763 255.1475 12.000000 55.32691 55.32691 2511.235 111.9716 114.67432 526.7413 548.4580 76.76997
Kumejima 597 -1.4178822 0.9250749 60.71743 54.01407 222.0000 277.0000 163.00000 2303.000 315.0000 132.5000 29.00000 809.5000 423.5000 576.5000 423.5000 52.00000 28.00000 4482.000 309.0000 128.50000 180.5000 260.0000 163.00000 1022.8400 0.5550000 12.77500 432.4150 95931.00 95931.00 182.0000 257.0000 12.000000 54.14000 54.14000 2539.940 111.3300 111.44500 532.5000 623.7271 119.42800
Miyakojima 552 -2.3345104 0.9870130 165.13666 42.67389 235.5000 281.2234 184.50000 2144.536 242.3613 136.8159 20.50000 646.8572 417.0039 621.6921 429.3743 49.00000 30.00000 3854.703 314.2234 153.89129 160.5653 270.2234 185.89129 1041.4816 0.5150000 10.99898 451.9025 101817.17 101817.17 201.1930 261.2234 12.000000 56.05820 66.00715 2539.382 111.8939 97.58236 590.6965 850.5334 120.19214
Ishigakijima 877 -0.6252511 0.7242757 230.48583 40.31216 236.1942 281.6206 185.25599 2132.099 236.0820 122.4948 21.94050 655.8379 388.7829 619.1044 400.0119 50.50000 31.00000 3833.381 313.2880 153.19419 160.0733 269.2880 187.28797 1055.7542 0.5071154 10.98712 447.8451 102097.75 102097.75 203.7891 263.2880 12.000000 58.52420 67.71221 2336.965 113.6409 99.85277 593.6196 954.2250 41.09759
Iriomotejima 887 0.1921938 0.3906094 292.10221 159.56565 228.9828 273.2398 178.98280 2370.300 258.4808 151.5683 18.23086 729.2650 471.4232 636.3513 479.9272 51.50000 32.12202 3736.810 305.5627 147.80150 158.2256 261.2398 183.77527 1050.9475 0.5553407 10.68895 506.5259 99032.11 99032.11 197.3049 252.8150 12.000000 57.64614 67.09907 2405.794 114.4487 99.20065 574.2911 989.3600 41.09759

12 Προκαταρτικές αναλύσεις

12.1 Έλεγχος κανονικότητας

12.1.1 Α

Δημιουργούμε ένα αρχείο το οποίο ελέγχει εάν όλες οι μεταβλητές μας (εκτός της πρώτης στήλης, η οποία είναι ποσοτική μεταβλητή - το όνομα των νησιών), έχουν κανονική κατανομή:

12.1.2 Β

Ας δούμε τα αποτελέσματα του προηγούμενου βήματος

## $SR
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.93472, p-value = 0.2893
## 
## 
## $PD
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.96192, p-value = 0.6968
## 
## 
## $PD_P
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.77558, p-value = 0.001314
## 
## 
## $Area
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.75456, p-value = 0.0007212
## 
## 
## $Elevation
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.72274, p-value = 0.0003039
## 
## 
## $bio1
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.93352, p-value = 0.2771
## 
## 
## $bio10
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.85281, p-value = 0.01499
## 
## 
## $bio11
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.93575, p-value = 0.3001
## 
## 
## $bio12
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.90384, p-value = 0.09261
## 
## 
## $bio13
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.86312, p-value = 0.02138
## 
## 
## $bio14
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92022, p-value = 0.1701
## 
## 
## $bio15
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92241, p-value = 0.1845
## 
## 
## $bio16
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.91524, p-value = 0.1414
## 
## 
## $bio17
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.88965, p-value = 0.055
## 
## 
## $bio18
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92892, p-value = 0.2344
## 
## 
## $bio19
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.89273, p-value = 0.06153
## 
## 
## $bio2
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.94329, p-value = 0.3913
## 
## 
## $bio3
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.91761, p-value = 0.1544
## 
## 
## $bio4
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.94292, p-value = 0.3863
## 
## 
## $bio5
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.82406, p-value = 0.005789
## 
## 
## $bio6
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.93515, p-value = 0.2938
## 
## 
## $bio7
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.94416, p-value = 0.4031
## 
## 
## $bio8
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.88975, p-value = 0.0552
## 
## 
## $bio9
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92902, p-value = 0.2352
## 
## 
## $annualPET
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.88922, p-value = 0.05414
## 
## 
## $climaticMoistureIndex
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92624, p-value = 0.2124
## 
## 
## $continentality
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.93748, p-value = 0.3193
## 
## 
## $embergerQ
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.91509, p-value = 0.1406
## 
## 
## $growingDegDays0
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92894, p-value = 0.2345
## 
## 
## $growingDegDays5
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92894, p-value = 0.2345
## 
## 
## $maxTempColdest
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92989, p-value = 0.2428
## 
## 
## $minTempWarmest
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.88613, p-value = 0.04841
## 
## 
## $monthCountByTemp10
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.46956, p-value = 1.199e-06
## 
## 
## $PETColdestQuarter
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9442, p-value = 0.4037
## 
## 
## $PETDriestQuarter
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92498, p-value = 0.2028
## 
## 
## $PETseasonality
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92548, p-value = 0.2065
## 
## 
## $PETWarmestQuarter
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.87995, p-value = 0.03873
## 
## 
## $PETWettestQuarter
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.91681, p-value = 0.1498
## 
## 
## $thermicityIndex
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.93209, p-value = 0.2631
## 
## 
## $MaxD
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.91265, p-value = 0.1284
## 
## 
## $MinD
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.8222, p-value = 0.005454

Οι μεταβλητές με πιθανότητα μικρότερη από 0.05 ΔΕΝ εμφανίζουν κανονική κατανομή

12.2 Έλεγχος συσχέτισης

Ας ελέγξουμε τη συσχέτιση μεταξύ των μεταβλητών μας
Εάν κάποιες έχουν συσχέτιση \(\geq 0.7\), πρέπει να τις αφαιρέσουμε από τη μεταγενέστερη ανάλυση (Γιατί;)

## Call:corr.test(x = ryukyus_new[, 5:42], method = "spearman")
## Correlation matrix 
##                        Area Elevation  bio1 bio10 bio11 bio12 bio13 bio14 bio15
## Area                   1.00      0.06  0.17  0.11  0.14 -0.01 -0.22  0.01 -0.32
##                       bio16 bio17 bio18 bio19  bio2  bio3  bio4  bio5  bio6
## Area                  -0.15 -0.07  0.25 -0.05 -0.22  0.11 -0.22  0.10  0.10
##                        bio7  bio8  bio9 annualPET climaticMoistureIndex
## Area                  -0.20  0.24  0.00      0.11                 -0.01
##                       continentality embergerQ growingDegDays0 growingDegDays5
## Area                           -0.22      0.32            0.16            0.16
##                       maxTempColdest minTempWarmest monthCountByTemp10
## Area                            0.15           0.09              -0.09
##                       PETColdestQuarter PETDriestQuarter PETseasonality
## Area                               0.18            -0.05          -0.44
##                       PETWarmestQuarter PETWettestQuarter thermicityIndex  MaxD
## Area                              -0.06             -0.23            0.14  0.04
##                        MinD
## Area                   0.49
##  [ reached getOption("max.print") -- omitted 37 rows ]
## Sample Size 
## [1] 16
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##                       Area Elevation bio1 bio10 bio11 bio12 bio13 bio14 bio15
## Area                  0.00      1.00 1.00  1.00  1.00  1.00  1.00  1.00  1.00
##                       bio16 bio17 bio18 bio19 bio2 bio3 bio4 bio5 bio6 bio7
## Area                   1.00  1.00  1.00  1.00 1.00 1.00 1.00 1.00 1.00 1.00
##                       bio8 bio9 annualPET climaticMoistureIndex continentality
## Area                  1.00 1.00      1.00                  1.00           1.00
##                       embergerQ growingDegDays0 growingDegDays5 maxTempColdest
## Area                       1.00            1.00            1.00           1.00
##                       minTempWarmest monthCountByTemp10 PETColdestQuarter
## Area                            1.00               1.00              1.00
##                       PETDriestQuarter PETseasonality PETWarmestQuarter
## Area                              1.00           1.00              1.00
##                       PETWettestQuarter thermicityIndex MaxD MinD
## Area                               1.00            1.00 1.00 1.00
##  [ reached getOption("max.print") -- omitted 37 rows ]
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

Διακρίνετε κάποιες μεταβλητές με τόσο μεγάλη συσχέτιση;

12.3 Οπτικοποίηση αποτελεσμάτων

Μπορούμε να οπτικοποιήσουμε τα αποτελέσματα μας και να διαπιστώσουμε εάν κάποια μεταβλητή δεν εμφανίζει κανονική κατανομή

12.4 Κανονικοποίηση των δεδομένων μας

Καθώς τα δεδομένα μας δεν έχουν κανονική κατανομή, θα πρέπει να τα κανονικοποιήσουμε. Αυτό σημαίνει ότι θα πρέπει να τα λογαριθμήσουμε.

13 Κυρίως ανάλυση

13.1 Μοντελοποίηση - ISAR

Τώρα είμαστε πλέον σε θέση να δημιουργήσουμε 2 μοντέλα απλής γραμμικής παλινδρόμησης, ένα για την ταξινομική και ένα για την φυλογενετική ποικιλότητα σε σχέση με την έκταση:

13.1.1 Περίληψη μοντέλων

Ας δούμε την περίληψη για κάθε ένα εκ των μοντέλων μας

## 
## Call:
## lm(formula = SR ~ Area, data = ryukyus_log)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.169622 -0.049189 -0.001934  0.072756  0.107531 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.17927    0.08117  26.848 1.93e-13 ***
## Area         0.28263    0.03750   7.537 2.72e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08492 on 14 degrees of freedom
## Multiple R-squared:  0.8023, Adjusted R-squared:  0.7882 
## F-statistic: 56.81 on 1 and 14 DF,  p-value: 2.722e-06
## 
## Call:
## lm(formula = PD ~ Area, data = ryukyus_log)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.02464 -0.50739  0.04005  0.60604  1.37789 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -2.4217     0.6971  -3.474  0.00372 **
## Area          0.5013     0.3221   1.557  0.14188   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7293 on 14 degrees of freedom
## Multiple R-squared:  0.1475, Adjusted R-squared:  0.08665 
## F-statistic: 2.423 on 1 and 14 DF,  p-value: 0.1419

ΠΡΟΣΟΧΗ
Μας ενδιαφέρει το adjusted R squared

13.1.2 Δημιουργία standardized residuals

Δημιουργούμε 2 αρχεία τα οποία περιέχουν τα standardized residuals των μοντέλων που δημιουργήσαμε στο προηγούμενο βήμα:

Μπορεί να γίνει σε ένα βήμα αντί για τρια:

## ============================================================================##
## Check the estimate and the p-value
## ============================================================================##
metrics <- lm(SR ~ Area, data = ryukyus_log) %>% tidy() %>% full_join(., lm(PD ~ 
    Area, data = ryukyus_log) %>% tidy()) %>% mutate(Metric = rep(c("SR", "PD"), 
    each = 2)) %>% dplyr::select(Metric, everything())
## ============================================================================##

## ============================================================================##
## Check the metrics
## ============================================================================##
metrics <- lm(SR ~ Area, data = ryukyus_log) %>% glance() %>% full_join(., lm(PD ~ 
    Area, data = ryukyus_log) %>% glance()) %>% mutate(Metric = c("SR", "PD")) %>% 
    dplyr::select(Metric, everything())
## ============================================================================##


## ============================================================================##
## Check the residuals
## ============================================================================##
resids <- lm(SR ~ Area, data = ryukyus_log) %>% augment() %>% mutate(Island = ryukyus_log$Island) %>% 
    rename(SR_res = .std.resid) %>% full_join(., lm(PD ~ Area, data = ryukyus_log) %>% 
    augment() %>% mutate(Island = ryukyus_log$Island) %>% rename(PD_res = .std.resid), 
    by = "Island") %>% dplyr::select(Island, SR, SR_res, PD, PD_res)
## ============================================================================##

13.1.3 Στικτόγραμμα

Τώρα θα δημιουργήσουμε ένα στικτόγραμμα για να δούμε ποιό νησί αποτελεί θερμό σημείο φυτικής ποικιλότητας για τον συνολικό αριθμό των φυτικών taxa, καθώς και για τον αριθμό των ενδημικών taxa:

13.2 Μοντελοποίηση - Πολλαπλή γραμμική παλινδρόμηση

13.2.1 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο enter

Τώρα είμαστε πλέον σε θέση να δημιουργήσουμε 2 μοντέλα πολλαπλής γραμμικής παλινδρόμησης, ένα για την ταξινομική και ένα για την φυλογενετική ποικιλότητα σε σχέση με όλες τις ανεξάρτητες, μη συγγραμικές μεταβλητές μας. Θα χρησιμοποιήσουμε αρχικά την μέθοδο enter:

ΠΡΟΣΟΧΗ
Τα ανωτέρω μοντέλα δεν δείχνουν ποιές είναι οι μεταβλητές οι οποίες επηρεάζουν τα πρότυπα ποικιλότητας. Δείχνουν απλά τι ποσοστό εξηγείται εάν συμπεριλαμβάνονται όλες οι ανεξάρτητες μεταβλητές στο μοντέλο μας.

Ας δούμε τώρα την περίληψη για κάθε ένα εξ αυτών των μοντέλων:

## 
## Call:
## lm(formula = SR ~ Area + bio14 + bio8 + monthCountByTemp10 + 
##     MaxD + MinD, data = ryukyus_log)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.135024 -0.049051  0.003223  0.031492  0.114373 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         1.09862    2.01050   0.546   0.5980    
## Area                0.27193    0.04050   6.715  8.7e-05 ***
## bio14               0.17172    0.48562   0.354   0.7318    
## bio8                1.04033    0.59001   1.763   0.1117    
## monthCountByTemp10 -1.52961    0.81899  -1.868   0.0946 .  
## MaxD                0.02294    0.30523   0.075   0.9417    
## MinD               -0.09494    0.13400  -0.708   0.4966    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08294 on 9 degrees of freedom
## Multiple R-squared:  0.8788, Adjusted R-squared:  0.7979 
## F-statistic: 10.87 on 6 and 9 DF,  p-value: 0.001092
## 
## Call:
## lm(formula = PD ~ Area + bio14 + bio8 + monthCountByTemp10 + 
##     MaxD + MinD, data = ryukyus_log)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8957 -0.2648 -0.0292  0.4049  0.7872 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)  
## (Intercept)        -20.5226    14.1621  -1.449   0.1812  
## Area                 0.7156     0.2853   2.509   0.0334 *
## bio14                7.5665     3.4207   2.212   0.0543 .
## bio8                 5.3023     4.1561   1.276   0.2340  
## monthCountByTemp10  -5.9614     5.7691  -1.033   0.3284  
## MaxD                -0.3244     2.1501  -0.151   0.8834  
## MinD                -2.1573     0.9439  -2.285   0.0481 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5842 on 9 degrees of freedom
## Multiple R-squared:  0.6484, Adjusted R-squared:  0.4139 
## F-statistic: 2.766 on 6 and 9 DF,  p-value: 0.08269

13.2.2 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο της βηματικής πολλαπλής γραμμικής παλινδρόμησης

Με την μέθοδο αυτή θα μπορέσουμε να δούμε ποιές μεταβλητές εμπερικλείονται στο βέλτιστο μοντέλο για τον συνολικό αριθμό των ειδών

Κάντε το ίδιο και για την φυλογενετική ποικιλότητα. Ποιές μεταβλητές εξηγούν και σε τι ποσοστό τα πρότυπα της φυλογενετικής ποικιλότητας;

ΠΡΟΣΟΧΗ: Μας ενδιαφέρει το adjusted R squared

## Start:  AIC=-12.41
## PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## - MaxD                1   0.00777 3.0795 -14.3653
## - monthCountByTemp10  1   0.36444 3.4361 -12.6119
## <none>                            3.0717 -12.4058
## - bio8                1   0.55550 3.6272 -11.7461
## - bio14               1   1.66988 4.7416  -7.4595
## - MinD                1   1.78274 4.8544  -7.0831
## - Area                1   2.14766 5.2194  -5.9234
## 
## Step:  AIC=-14.37
## PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## - monthCountByTemp10  1   0.38000 3.4595 -14.5036
## <none>                            3.0795 -14.3653
## - bio8                1   0.56038 3.6398 -13.6903
## + MaxD                1   0.00777 3.0717 -12.4058
## - bio14               1   1.73902 4.8185  -9.2021
## - MinD                1   1.85944 4.9389  -8.8071
## - Area                1   2.13989 5.2194  -7.9234
## 
## Step:  AIC=-14.5
## PD ~ Area + bio14 + bio8 + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## - bio8                1   0.33274 3.7922 -15.0342
## <none>                            3.4595 -14.5036
## + monthCountByTemp10  1   0.38000 3.0795 -14.3653
## + MaxD                1   0.02333 3.4361 -12.6119
## - bio14               1   1.46668 4.9261 -10.8485
## - MinD                1   1.81014 5.2696  -9.7701
## - Area                1   2.38991 5.8494  -8.1001
## 
## Step:  AIC=-15.03
## PD ~ Area + bio14 + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## <none>                            3.7922 -15.0342
## + bio8                1   0.33274 3.4595 -14.5036
## + monthCountByTemp10  1   0.15236 3.6398 -13.6903
## + MaxD                1   0.00319 3.7890 -13.0477
## - bio14               1   1.21835 5.0106 -12.5767
## - MinD                1   1.62428 5.4165 -11.3303
## - Area                1   2.38796 6.1802  -9.2199
## 
## Call:
## lm(formula = PD ~ Area + bio14 + MinD, data = ryukyus_log)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.99647 -0.28430 -0.05602  0.32450  0.81760 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -12.8387     6.5774  -1.952   0.0747 .
## Area          0.7480     0.2721   2.749   0.0176 *
## bio14         5.8728     2.9910   1.963   0.0732 .
## MinD         -1.5102     0.6661  -2.267   0.0427 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5622 on 12 degrees of freedom
## Multiple R-squared:  0.5659, Adjusted R-squared:  0.4573 
## F-statistic: 5.214 on 3 and 12 DF,  p-value: 0.01554
## Start:  AIC=-12.41
## PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## - MaxD                1   0.00777 3.0795 -14.3653
## - monthCountByTemp10  1   0.36444 3.4361 -12.6119
## <none>                            3.0717 -12.4058
## - bio8                1   0.55550 3.6272 -11.7461
## - bio14               1   1.66988 4.7416  -7.4595
## - MinD                1   1.78274 4.8544  -7.0831
## - Area                1   2.14766 5.2194  -5.9234
## 
## Step:  AIC=-14.37
## PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## - monthCountByTemp10  1   0.38000 3.4595 -14.5036
## <none>                            3.0795 -14.3653
## - bio8                1   0.56038 3.6398 -13.6903
## + MaxD                1   0.00777 3.0717 -12.4058
## - bio14               1   1.73902 4.8185  -9.2021
## - MinD                1   1.85944 4.9389  -8.8071
## - Area                1   2.13989 5.2194  -7.9234
## 
## Step:  AIC=-14.5
## PD ~ Area + bio14 + bio8 + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## - bio8                1   0.33274 3.7922 -15.0342
## <none>                            3.4595 -14.5036
## + monthCountByTemp10  1   0.38000 3.0795 -14.3653
## + MaxD                1   0.02333 3.4361 -12.6119
## - bio14               1   1.46668 4.9261 -10.8485
## - MinD                1   1.81014 5.2696  -9.7701
## - Area                1   2.38991 5.8494  -8.1001
## 
## Step:  AIC=-15.03
## PD ~ Area + bio14 + MinD
## 
##                      Df Sum of Sq    RSS      AIC
## <none>                            3.7922 -15.0342
## + bio8                1   0.33274 3.4595 -14.5036
## + monthCountByTemp10  1   0.15236 3.6398 -13.6903
## + MaxD                1   0.00319 3.7890 -13.0477
## - bio14               1   1.21835 5.0106 -12.5767
## - MinD                1   1.62428 5.4165 -11.3303
## - Area                1   2.38796 6.1802  -9.2199
## 
## Call:
## lm(formula = PD ~ Area + bio14 + MinD, data = ryukyus_log)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.99647 -0.28430 -0.05602  0.32450  0.81760 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -12.8387     6.5774  -1.952   0.0747 .
## Area          0.7480     0.2721   2.749   0.0176 *
## bio14         5.8728     2.9910   1.963   0.0732 .
## MinD         -1.5102     0.6661  -2.267   0.0427 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5622 on 12 degrees of freedom
## Multiple R-squared:  0.5659, Adjusted R-squared:  0.4573 
## F-statistic: 5.214 on 3 and 12 DF,  p-value: 0.01554

13.2.2.1 Ποιό μοντέλο είναι καλύτερο;

Με την βοήθεια του Akaike Information Criterion (εν συντομία: AIC), μπορούμε να διαπιστώσουμε ποιό μοντέλο είναι “καλύτερο”. Όσο μικρότερη η τιμή του AIC, τόσο καλύτερο το μοντέλο.

Κάντε το ίδιο και για την φυλογενετική ποικιλότητα. Ποιό από τα δύο μοντέλα για κάθε κατηγορία ειδών είναι καλύτερο;

13.2.2.3 Υπολογισμός της σχετικής σημαντικότητας των μεταβλητών

Μετά την περισπωμένη (~) βάζουμε όποιες μεταβλητές έχουν βγει σημαντικές από την βηματική πολλαπλή γραμμική παλινδρόμηση:

Weights
MaxD 4.671829
bio8 4.868050
monthCountByTemp10 5.001657
MinD 22.456880
bio14 31.196091
Area 31.805492

13.2.3 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο των ολικών υποσυνόλων

Μέχρι στιγμής, έχουμε δει δύο διαφορετικές μεθόδους, όσον αφορά την πολλαπλή γραμμική παλινδρόμηση. Όμως, δεν μπορούμε να είμαστε 100% σίγουροι ότι έχουμε φθάσει στο ορθό συμπέρασμα, καθότι ένα από τα μειονεκτήματα της βηματικής πολλαπλής παλινδρόμησης, είναι ότι μπορεί να “κολλήσει” σε ένα υποβέλτιστο πλατώ. Τον σκόπελο αυτό, μπορούμε να τον αποφύγουμε χρησιμοποιώντας την μέθοδο των ολικών υποσυνόλων:

13.2.3.1 Απεικόνιση των αποτελεσμάτων του προηγούμενου βηματος

Το γράφημα αυτό θα μας δείξει ποιός είναι ο πραγματικά βέλτιστος συνδυασμός των ανεξάρτητων μεταβλητών:

## 
## Call:
## lm(formula = PD ~ Area + bio14 + bio8 + monthCountByTemp10 + 
##     MinD, data = ryukyus_log)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.89706 -0.24475 -0.02862  0.37684  0.80405 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)  
## (Intercept)        -21.6286    11.5101  -1.879   0.0897 .
## Area                 0.7130     0.2705   2.636   0.0249 *
## bio14                7.4248     3.1244   2.376   0.0389 *
## bio8                 5.3227     3.9457   1.349   0.2071  
## monthCountByTemp10  -5.6306     5.0687  -1.111   0.2926  
## MinD                -2.1181     0.8620  -2.457   0.0338 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5549 on 10 degrees of freedom
## Multiple R-squared:  0.6475, Adjusted R-squared:  0.4712 
## F-statistic: 3.673 on 5 and 10 DF,  p-value: 0.03798
Weights
bio8 5.224379
monthCountByTemp10 5.706259
MinD 24.302432
Area 31.886016
bio14 32.880914

13.2.4 Έλεγχος των αποτελεσμάτων της παλινδρόμησης

Τώρα θα ελέγξουμε κατ’αρχάς εάν οι μεταβλητές του βέλτιστου μοντέλου μας παραβιάζουν μια από τις αρχές της πολλαπλής γραμμικής παλινδρόμησης:

##               Area              bio14               bio8 monthCountByTemp10 
##              FALSE              FALSE              FALSE              FALSE 
##               MinD 
##              FALSE

13.2.5 Ποιό μοντέλο είναι καλύτερο;

Με την βοήθεια του Akaike Information Criterion (εν συντομία: AIC), μπορούμε να διαπιστώσουμε ποιό μοντέλο είναι “καλύτερο”. Όσο μικρότερη η τιμή του AIC, τόσο καλύτερο το μοντέλο.

df AIC
enter2 8 35.00027
stepm2 5 32.37180
tot 7 33.04070

13.2.6 Έλεγχος πιθανοτήτων

Στο σημείο αυτό θα χρειαστεί να ελέγξουμε εάν η πιθανότητα κάθε μεταβλητής η οποία εμπερικλείεται στο τελικό μας μοντέλο είναι πραγματικά στατιστικώς σημαντική:

## 
## Call:
## lm(formula = PD ~ Area + bio14 + bio8 + monthCountByTemp10 + 
##     MinD, data = ryukyus_log)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.89706 -0.24475 -0.02862  0.37684  0.80405 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)  
## (Intercept)        -21.6286    11.5101  -1.879   0.0897 .
## Area                 0.7130     0.2705   2.636   0.0249 *
## bio14                7.4248     3.1244   2.376   0.0389 *
## bio8                 5.3227     3.9457   1.349   0.2071  
## monthCountByTemp10  -5.6306     5.0687  -1.111   0.2926  
## MinD                -2.1181     0.8620  -2.457   0.0338 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5549 on 10 degrees of freedom
## Multiple R-squared:  0.6475, Adjusted R-squared:  0.4712 
## F-statistic: 3.673 on 5 and 10 DF,  p-value: 0.03798
## [1] 0.1245 0.1945 1.0000 1.0000 0.1690

14 Θερμά σημεία βιοποικιλότητας

14.1 Εισαγωγή

Οι κλιμακούμενες απειλές για την βιοποικιλότητα, σε συνδυασμό με την υποχρηματοδότηση, έχουν ως αποτέλεσμα την εφαρμογή μιας διαδικασίας προτεραιοποίησης όσον αφορά τη λήψη αποφάσεων διαχειριστικού χαρακτήρα. Η διαδικασία αυτή επικεντρώνεται στην ιεράρχηση των ειδών, των πληθυσμών ή/και των οικοτόπων, με βάση τα οφέλη της βιοποικιλότητας, την δυνατότητα ανάκαμψης των ειδών και το κόστος επίτευξης του επιθυμητού στόχου (Farrier et al., 2007; Bottrill et al., 2008; Arponen, 2012; Jetz et al., 2014; Reece & Noss, 2014).

Ένα αποτελεσματικό εργαλείο για αυτόν τον σκοπό, είναι ο εντοπισμός θερμών σημείων ταξινομικής και φυλογενετικής βιοποικιλότητας.

Η κύρια στρατηγική για τη διατήρηση της βιοποικιλότητας έγκειται στον προσδιορισμό περιοχών για την προάσπιση των οργανισμών έναντι της εκάστοτε υφιστάμενης ή μελλοντικής απειλής (Margules & Pressey, 2000). Τα δίκτυα προστατευόμενων περιοχών (ΠΠ) αποτελούν τη ναυαρχίδα των πολιτικών διατήρησης της βιοποικιλότητας, ενώ διαδραματίζουν σημαίνοντα ρόλο στην προστασία της (Chape et al., 2005), καθώς η διατήρηση αρκετών ειδών είναι ισχυρώς εξαρτώμενη από τα δίκτυα αυτά (Watson et al., 2014). Παγκοσμίως, σχεδόν το 13% της χερσαίας έκτασης βρίσκεται υπό κάποιο καθεστώς προστασίας, ενώ το αντίστοιχο ποσοστό στις χώρες του ανεπτυγμένου κόσμου ανέρχεται στο ~12%, με εμφανή τάση αύξησης τα τελευταία χρόνια (Jenkins & Joppa, 2009). Για την αντιμετώπιση της παρατηρούμενης τάσης μείωσης της βιοποικιλότητας, η Συνθήκη για τη Βιολογική Ποικιλότητα, μέσω της Παγκόσμιας Στρατηγικής για τη Διατήρηση της Φυτικής και Ζωικής Ποικιλότητας, έχει θέσει ως στόχο την ανάσχεση της μείωσης αυτής (Jackson & Kennedy, 2009), δια της επίτευξης των Aichi Targets.

Η προτεραιοποίηση των περιοχών διατήρησης εδράζεται στον εντοπισμό θερμών σημείων βιοποικιλότητας (ΘΣΒ), τα οποία παραδοσιακά προσδιορίζονται βάσει ταξινομικών δεικτών [π.χ., αριθμός ειδών/ενδημικών, ποσοστό ενδημισμού, πιθανότητα εξαφάνισης – Forest et al. (2007)] και αποτελούν σημαντικά εργαλεία στην διαδικασία διαλογής όσον αφορά τη λήψη αποφάσεων διαχειριστικού χαρακτήρα. Ωστόσο, η προσέγγιση αυτή δεν λαμβάνει υπόψη ότι όλα τα είδη δεν είναι ισοδύναμα όσον αφορά την εξελικτική τους ιστορία ή/και τις λειτουργίες που επιτελούν (Cadotte & Jonathan Davies, 2010), ούτε το γεγονός ότι οι διαειδικές φυλογενετικές σχέσεις και κατ’ επέκταση η φυλογενετική ποικιλότητα της εκάστοτε περιοχής, αντικατοπτρίζουν τις οικολογικές, εξελικτικές και βιογεωγραφικές διαδικασίες με τις οποίες παράγεται, κατανέμεται και διατηρείται η βιοποικιλότητα (Cadotte & Jonathan Davies, 2010; Daru et al., 2019). Δεδομένου ότι ορισμένα είδη έχουν περισσότερο διακριτή εξελικτική ιστορία, η μη τυχαία εξαφάνιση ειδών μπορεί να οδηγήσει ορισμένους φυλογενετικούς κλάδους να απωλέσουν μεγαλύτερο ποσοστό ειδών σε σχέση με άλλους (Davies & Yessoufou, 2013). Καθώς είναι πλέον εφικτή η ενσωμάτωση φυλογενετικής πληροφορίας (Smith & Brown, 2018) στον προσδιορισμό και την προτεραιοποίηση των περιοχών διατήρησης, το γεγονός αυτό μπορεί να αποτελέσει ένα από τα σημαντικότερα βήματα προς την αποτελεσματικότερη διατήρηση και προστασία της βιοποικιλότητας (Tucker et al., 2019) σε διάφορες χωρικές και ταξινομικές κλίμακες (Pollock et al., 2017).

Στο φυλογενετικό επίπεδο, οι ισοδύναμοι των παραδοσιακών, ταξινομικών δεικτών για τον εντοπισμό ΘΣΒ, είναι η φυλογενετική ποικιλότητα [ΦΠ, sensu Faith (1992) – ισοδύναμη του αριθμού των ειδών], ο φυλογενετικός ενδημισμός [ΦΕ, sensu Rosauer et al. (2009) – ισοδύναμος του αριθμού των ενδημικών/ποσοστού ενδημισμού] και ο δείκτης EDGE [Evolutionary Distinct and Globally Endangered – EDGE, sensu Isaac et al. (2007) – ισοδύναμος της πιθανότητας εξαφάνισης]. Οι εν λόγω δείκτες ποσοτικοποιούν διαφορετικές πτυχές της εξελικτικής ποικιλότητας (Tucker et al., 2017).

Η ΦΠ είναι το άθροισμα του μήκους των κλάδων που συνδέουν ένα σύμπλεγμα ειδών με την ρίζα ενός φυλογενετικού δένδρου (Faith, 1992).

Ο ΦΕ ποσοτικοποιεί τον βαθμό στον οποίο ένα σημαντικό ποσοστό της ΦΠ εντοπίζεται αποκλειστικά εντός της εκάστοτε περιοχής μελέτης (Rosauer et al., 2009).

Ο δείκτης EDGE συνδυάζει την εξελικτική διακριτότητα (ED, η φυλογενετική απομόνωση ενός είδους) με το καθεστώς κινδύνου εξαφάνισης σε παγκόσμιο επίπεδο (GE) σύμφωνα με τα πρότυπα της IUCN (Isaac et al., 2007).

Οι δείκτες ΦΕ και EDGE αντιπροσωπεύουν σταθμισμένες παραλλαγές της ΦΠ σε γεωγραφική κλίμακα και καθεστώς απειλής, αντίστοιχα (Daru et al., 2019). Περιοχές με στατιστικώς σημαντικά υψηλό ή χαμηλό ΦΕ είναι κέντρα πάλαιο- ή νέο-ενδημισμού, αντίστοιχα (Mishler et al., 2014).

Ένας ακόμα σημαντικός δείκτης για τον ασφαλή εντοπισμό ΘΣΒ, είναι η σταθμισμένη ταξινομική ποικιλότητα (Corrected weighted endemism), καθώς δίνει περισσότερη έμφαση στα στενότοπα taxa (Crisp et al., 2001; Linder, 2001a, 2001b)

Συνεπώς, συνδυάζοντας τους ταξινομικούς και φυλογενετικούς αυτούς δείκτες, ο εντοπισμός των ΘΣΒ και κατ’ επέκταση η προτεραιοποίηση των ΠΠ αποκτά έναν πιο ολοκληρωμένο και ολιστικό χαρακτήρα που πιθανότατα επαρκεί για τη μακροπρόθεσμη διατήρηση της βιοποικιλότητας (Buerki et al., 2015), μειώνοντας το οικονομικό κόστος των μέτρων διατήρησης (Carta et al., 2019).

14.2 Εντοπισμός θερμών σημείων βιοποικιλότητας

Σε αυτό το σημείο θα χρησιμοποιήσουμε τα δεδομένα από την Ιταλία, τα οποία δημιουργήσαμε σε προηγούμενο εργαστήριο.

14.2.1 Εισαγωγή δεδομένων

Ας εισαγάγουμε τα δεδομένα για την Ιταλία.

Θυμηθείτε ότι θα χρειαστεί να αλλάξετε το file path και να χρησιμοποιήσετε τα δεδομένα από το Github.

14.2.4 Δημιουργία μήτρας παρουσίας/απουσίας

Και ακολούθως ένα αρχείο το οποίο θα αποτελέσει την βάση των αναλύσεων μας:

Ας δούμε πως είναι:

## 5 x 5 sparse Matrix of class "dgCMatrix"
##        Abies_alba Abies_nebrodensis Abies_nordmanniana Abutilon_theophrasti
## v10000          .                 .                  .                    .
## v10003          .                 .                  .                    .
## v10004          .                 .                  .                    .
## v10007          .                 .                  .                    .
## v10008          .                 .                  .                    .
##        Acacia_dealbata
## v10000               .
## v10003               .
## v10004               .
## v10007               .
## v10008               .

Και ας το αποθηκεύσουμε:

14.2.5 Κατασκευή φυλογενετικού δένδρου

Ας κατασκευάσουμε το φυλογενετικό δένδρο για τα taxa τα οποία απαντώνται στην Ιταλία.

Υπάρχει πιθανότητα αναλόγως των δυνατοτήτων του Η/Υ σας, η εντολή αυτή να διαρκέσει μεγάλο χρονικό διάστημα. Μπορείτε να εισαγάγετε το αποτέλεσμα αυτής με την εντολή: italy_tree <- readRDS(Italian phylogenetic tree.rds.

14.2.6 Κυρίως ανάλυση

Τώρα μπορούμε να δημιουργήσουμε ένα αρχείο το οποίο θα περιέχει πληροφορία σχετικά με το ποια κελιά αποτελούν θερμά σημεία σταθμισμένης και μη-σταθμισμένης ταξινομικής και φυλογενετικής ποικιλότητας:

14.2.7 Οπτικοποίηση

Ας δημιουργήσουμε ένα νέο αρχείο το οποίο θα χρησιμοποιήσουμε για την οπτικοποίηση των αποτελέσματων μας και το οποίο θα περιέχει μια νέα μεταβλητή, η οποία θα υποδεικνύει ποια κελιά αποτελούν θερμά ή ψυχρά σημεία βιοποικιλότητας για όλους τους δείκτες που χρησιμοποιήσαμε:

Απεικονίστε τα θερμά σημεία βιοποικιλότητας σύμφωνα με τους δείκτες PD & PE.

15 Appendix: R code

##===========================================================================##
## Do NOT run this code
##===========================================================================##
## Global options
options(Encoding="UTF-8")

library(knitr)
library(rmdformats)

## Global options
options(max.print="75")
opts_chunk$set(echo=TRUE,
                 cache=TRUE,
               prompt=FALSE,
               tidy=TRUE,
               # comment=NA,
               message=FALSE,
               warning=FALSE,
               eval = TRUE)
opts_knit$set(width=75)
##===========================================================================##
htmltools::img(src = knitr::image_uri(file.path("G:/Academic/Lessons/EMB/Labs/Labs/TRUE RMD Files/HBS Seminar/Logo_HBS/ebe-LOGO-005-300_1_70.jpg")), 
               alt = 'logo', 
               style = 'position:absolute; top:0; right:0; padding:5px;')
##===========================================================================##
## Load the libraries
##===========================================================================##
library(tidyverse)
library(raster)
library(sf)
library(ape)
library(phyloregion)
library(picante)
library(rtrees)
library(ggtree)
library(exactextractr)
library(caret)
library(psych)
library(PhyloMeasures)
library(phyloregion)
library(Matrix)
library(broom)
##===========================================================================##
##===========================================================================##
## Load our data
##===========================================================================##
plants <- read.csv("G:/Academic/Lessons/EMB/Labs/Labs/Ryukyus.csv", h = T, sep = ";", row.names = 1)
##===========================================================================##
##===========================================================================##
## Check their dimensions
##===========================================================================##
dim(plants) 
##===========================================================================##
## Check their dimensions
##===========================================================================##
str(plants) 
pre {
  max-height: 400px;
  overflow-y: auto;
}

pre[class] {
  max-height: 400px;
}
##===========================================================================##
## Convert the P/A matrix to SR
##===========================================================================##
ryukyus <- rowSums(plants) %>% 
  enframe() %>% 
  dplyr::rename(Island = name, 
                SR = value)

ryukyus
##===========================================================================##
## Create a vector containing just the file paths
##===========================================================================##
shps <- fs::dir_ls("C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Data/Vectors/GIS/JAPAN/", 
                   glob = "*.shp")
##===========================================================================##


##===========================================================================##
## Create a vector containing just the island names
##===========================================================================##
island_nms <- shps %>% 
  str_replace_all('C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Data/Vectors/GIS/JAPAN/',
                  '') %>% 
  str_replace_all('\\.shp', '')
##===========================================================================##


##===========================================================================##
## Bulk load the shapefiles
##===========================================================================##
ryukyus_shps <- shps %>% 
  purrr::map(read_sf) %>% 
  lapply(., dplyr::select, NAME_1, geometry)  %>% 
  do.call(rbind, .) %>% 
  mutate(Island = island_nms) %>% 
  dplyr::select(Island, geometry)
##===========================================================================##
##===========================================================================##
## Load, crop and mask the rasters
##===========================================================================##
preds <- stack(stack(list.files(path = "I:/SDM DATA/WORLDCLIM/Current/bio_2-5m_bil/",
                                pattern = '.bil',
                                recursive = T,
                                full.names = T)) %>% 
                 crop(., ryukyus_shps, snap = 'in') %>% 
                 mask(., ryukyus_shps),
               stack(list.files(path = "I:/SDM DATA/Envirem_Predictors/Eurasia_current_2.5arcmin_geotiff/",
                                pattern = '.tif',
                                recursive = T,
                                full.names = T))  %>% 
                 crop(., ryukyus_shps, snap = 'in') %>% 
                 mask(., ryukyus_shps))
##===========================================================================##
##===========================================================================##
## Load, crop and mask the altitudinal raster
##===========================================================================##
altitude <- raster("I:/SDM DATA/Altitude/elevation_1KMmd_GMTEDmd.tif")  %>% 
  crop(., ryukyus_shps, snap = 'in') %>% 
  mask(., ryukyus_shps)
##===========================================================================##
##===========================================================================##
## Extract abiotic info
##===========================================================================##
preds_ryukyus <- exact_extract(preds, ryukyus_shps, 'median') %>% 
  mutate(Island = ryukyus_shps$Island) %>% 
  full_join(.,
            exact_extract(altitude, ryukyus_shps, 'median') %>% 
              enframe() %>% 
              mutate(Island = ryukyus_shps$Island) %>% 
              rename(Elevation = value) %>% 
              dplyr::select(Island, Elevation)) %>% 
  dplyr::select(Island, everything()) %>% 
  as_tibble()


##============================================================================##
## Give them proper names
##============================================================================##
col_names <- preds_ryukyus %>%
  colnames() %>%
  enframe() %>%
  mutate(value = str_replace_all(value, 'median\\.', '')) %>%
  mutate(value = str_replace_all(value, 'current_2\\.5arcmin_', '')) %>% 
  pull(value)
##===========================================================================##


##============================================================================##
## Overwrite names and remove Row 1
##============================================================================##
vars <- preds_ryukyus %>%
  set_names(col_names) %>% 
  dplyr::select(Island, everything()) %>% 
  as_tibble()
##============================================================================##
## Transform the missing data
##============================================================================##
imputed_vars <- predict(preProcess(as.data.frame(vars[, 2:36]), 
                                   method = c('medianImpute')), 
                        vars)
##============================================================================##
## First drop any rows with NA values
##============================================================================##
drop_vars <- vars %>% drop_na()
##============================================================================##


##============================================================================##
## Then check which islands have NA values in any of the variables
##============================================================================##
setdiff(imputed_vars$Island, drop_vars$Island)


##============================================================================##
## Create a symbol to filter out all the Islands having NA values
##============================================================================##
'%!in%' <- function(x,y)!('%in%'(x,y))
##============================================================================##


##============================================================================##
## Now create a temporary presence/absence matrix removing any problematic islands
##============================================================================##
plants_new <- plants %>% 
  mutate(Island = rownames(.)) %>% 
  as_tibble() %>% 
  dplyr::select(Island, everything()) %>% 
  pivot_longer(!Island, names_to = 'Taxon', 
               values_to = 'PA') %>% 
  filter(Island %!in% setdiff(imputed_vars$Island, 
                              drop_vars$Island)) %>% 
  pivot_wider(names_from = Taxon,
              values_from = PA)
##============================================================================##


##============================================================================##
## See which taxa are left after removing the problematic islands
##============================================================================##
keepers <- colSums(plants_new[, -1]) %>% 
  enframe() %>% 
  filter(value > 0)
##============================================================================##


##============================================================================##
## Finally create the new PA matrix
##============================================================================##
plants_new_pa <- plants_new %>% 
  pivot_longer(!Island, names_to = 'Taxon', 
               values_to = 'PA') %>% 
  filter(Taxon %in% keepers$name) %>% 
  pivot_wider(names_from = Taxon,
              values_from = PA)
##============================================================================##
datatable(plants, options = list(pagelength = 10)) 
jpn_tree <- readRDS('G:/Academic/Lessons/EMB/Labs/Labs/16 Ryukyus tree.rds')
##============================================================================##
## Plot the tree
##============================================================================##
groupInfo <- split(jpn_tree$tip.label, 
                   gsub("_\\w+", "", 
                        jpn_tree$tip.label))

groups_jpn <- groupOTU(jpn_tree, 
                       groupInfo)
ggtree(groups_jpn, aes(color = group),
       layout = 'circular') +
  geom_tiplab2(size = 0.5, align = T) 
##============================================================================##
## Combine the trees with the P/A data
##============================================================================##

#check for mismatches/missing species
combined_japan <- picante::match.phylo.comm(jpn_tree,plants_new_pa[, 2:1800])

##============================================================================##
## The resulting object is a list with $phy and $comm elements
## replace our original data with the sorted/matched data
##============================================================================##
phy_japan <- combined_japan$phy

comm_japan <- combined_japan$comm
##============================================================================##
##============================================================================##
## Compute phylogenetic alpha diversity
##============================================================================##

## For the estimation of SES values
japan_pd <- pd.query(phy_japan, comm_japan, TRUE) 

## For the p-values
japan_pd_p <- pd.pvalues(phy_japan, comm_japan, reps = 1000)
##============================================================================##
##============================================================================##
## Create the final dataset
##============================================================================##
ryukyus_new <- ryukyus %>% 
  filter(Island %in% plants_new_pa$Island) %>% 
  full_join(., drop_vars) %>% 
  full_join(., ryukyus_shps) %>% 
  filter(Island %!in% setdiff(imputed_vars$Island, 
                              drop_vars$Island)) %>% 
  st_as_sf() %>% 
  
  ## Lambert Azimuthal Equal Area projection
  st_transform('+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs') %>% 
  
  ## Create the Area and PD variables
  mutate(Area = as.numeric(units::set_units(st_area(.), km^2)),
         PD = japan_pd,
         PD_P = japan_pd_p) %>% 
  
  ## Transform back to WGS84
  st_transform(4326) %>% 
  
  dplyr::select(Island, SR, PD, PD_P, Area, Elevation, everything()) %>% 
  
  ## It has zero-variance
  dplyr::select(-aridityIndexThornthwaite) %>% 
  
  ## Drop the geometry
  st_drop_geometry() %>% 
  
  ## Join these data with the Maximum and Minimum distance between each island
  ## pair
  inner_join(., 
             
             ## Calculate the Maximum distance between each island pair
             apply((st_distance(ryukyus_shps %>% 
                                     st_centroid(), 
                                   ryukyus_shps %>% 
                                     st_centroid())/1e+3) %>% 
                        as_tibble() %>% 
                        set_names(ryukyus_shps$Island) %>% 
                        units::drop_units() %>% 
                        na_if(., 0), 1, max, na.rm = T) %>% 
               enframe() %>% 
               rename(MaxD = value) %>% 
               mutate(Island = ryukyus_shps$Island) %>% 
               dplyr::select(Island, MaxD) %>% 
               
               full_join(., 
                         
                         ## Calculate the Minimum distance between each island
                         ## pair
                         apply((st_distance(ryukyus_shps %>% 
                                                 st_centroid(), 
                                               ryukyus_shps %>% 
                                                 st_centroid())/1e+3) %>% 
                                    as_tibble() %>% 
                                    set_names(ryukyus_shps$Island) %>% 
                                    units::drop_units() %>% 
                                    na_if(., 0), 1, min, na.rm = T) %>% 
                           enframe() %>% 
                           rename(MinD = value) %>% 
                           mutate(Island = ryukyus_shps$Island) %>% 
                           dplyr::select(Island, MinD)))
##============================================================================##

ryukyus_new
##============================================================================##
## Create the phylogenetic tree
##============================================================================##
library(rtrees)
jpn_tree <- keepers %>% 
  separate(name, into = c('genus', 'binomial'), sep = '_', remove = F) %>% 
  rename(species = name) %>% 
  dplyr::select(species, genus) %>% 
  get_tree(sp_list = .,
           tree = tree_plant_otl, # either 
           taxon = "plant", # or
           scenario = "S1",
           show_grafted = FALSE)
##============================================================================##
##============================================================================##
## Normality check
##============================================================================##
normality <- lapply(ryukyus_new[,2:42], shapiro.test)
##============================================================================##
normality
##============================================================================##
## Correlation check
##============================================================================##
corr.test(ryukyus_new[,5:42], method = "spearman")
##============================================================================##
## Visualization
##============================================================================##
car::scatterplotMatrix(ryukyus_new[,5:10], 
                  spread = F,
                  smoother.args = list(lty = 2), 
                  main = "Scatter Plot Matrix")
##============================================================================##
## Correlation check
##============================================================================##
unconc <- usdm::vifcor(ryukyus_new %>% 
                         as.data.frame() %>% 
                         dplyr::select(Area:MinD), 
                       th = 0.7)
unconc <- unconc@results$Variables[unconc@results$VIF <= 5]
unconc <- droplevels(factor(unconc))
##============================================================================##
##============================================================================##
## Log-transform the data
##============================================================================##
ryukyus_log <- ryukyus_new %>% 
  dplyr::select(Area, bio14, bio8,
                monthCountByTemp10, 
                MaxD, MinD) %>% 
  mutate(SR = ryukyus_new$SR) %>% 
  log10(.) %>% 
  mutate(Island = ryukyus_new$Island,
         PD = ryukyus_new$PD) %>% 
  dplyr::select(Island, SR, PD, everything())
##============================================================================##
##============================================================================##
## ISAR models
##============================================================================##
model1 <- lm(SR ~ Area, data = ryukyus_log)
model2 <- lm(PD ~ Area, data = ryukyus_log)
##============================================================================##
##============================================================================##
## Check the summary
##============================================================================##
summary(model1)
##============================================================================##
## Check the summary
##============================================================================##
summary(model2)
##============================================================================##
## Create the residuals
##============================================================================##
res_sr <- rstandard(model1)
res_pd <- rstandard(model2)
##============================================================================##
##============================================================================##
## Check the estimate and the p-value
##============================================================================##
metrics <- lm(SR ~ Area, data = ryukyus_log) %>% tidy() %>% 
  full_join(., lm(PD ~ Area, data = ryukyus_log) %>% tidy()) %>% 
  mutate(Metric = rep(c('SR','PD'), each = 2)) %>% 
  dplyr::select(Metric, everything())
##============================================================================##

##============================================================================##
## Check the metrics
##============================================================================##
metrics <- lm(SR ~ Area, data = ryukyus_log) %>% glance() %>% 
  full_join(., lm(PD ~ Area, data = ryukyus_log) %>% glance()) %>% 
  mutate(Metric = c('SR','PD')) %>% 
  dplyr::select(Metric, everything())
##============================================================================##


##============================================================================##
## Check the residuals
##============================================================================##
resids <- lm(SR ~ Area, 
             data = ryukyus_log) %>% 
  augment() %>% 
  mutate(Island = ryukyus_log$Island)  %>% 
  rename(SR_res = .std.resid) %>% 
  full_join(., lm(PD ~ Area, data = ryukyus_log) %>% 
              augment() %>% 
              mutate(Island = ryukyus_log$Island) %>% 
              rename(PD_res = .std.resid),
            by = 'Island') %>% 
  dplyr::select(Island, SR, SR_res, PD, PD_res)
##============================================================================##
##============================================================================##
## Plot
##============================================================================##
ggplot(ryukyus_log, aes(res_sr, res_pd)) +
  geom_point() + 
  geom_text(aes(label = Island), 
            fontface = "bold", 
            size = 2.8, 
            vjust = -1.2, 
            hjust = 0.5) + 
  geom_hline(yintercept = 0.0) + 
  geom_vline(xintercept = 0.0) + 
  theme(legend.position = "none") + 
  xlab("SR") + 
  ylab("PD") 
##============================================================================##
## Quick save
##============================================================================##
ggsave("SR-PD.png", width = 20, height = 20, units = "cm", dpi = 600)
##============================================================================##
## The enter method
##============================================================================##
enter1 <- lm(SR ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log) 

enter2 <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log) 
##============================================================================##
##============================================================================##
## Check the summary
##============================================================================##
summary(enter1)
summary(enter2)

##============================================================================##
## The solution
##============================================================================##
stepm2 <- MASS::stepAIC(lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log),
                        scale = 0,
                        direction = "both",
                        trace = 1,
                        keep = NULL,
                        steps = 1000,
                        k = 2)

summary(stepm2)
##============================================================================##
## The stepAIC method
##============================================================================##
stepm1 <- MASS::stepAIC(lm(SR ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log), 
                  scale = 0, 
                  direction = "both", 
                  trace = 1, 
                  keep = NULL, 
                  steps = 1000,
                  k = 2)

stepm2 <- MASS::stepAIC(lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log),
                        scale = 0,
                        direction = "both",
                        trace = 1,
                        keep = NULL,
                        steps = 1000,
                        k = 2)

summary(stepm1)
stepm2 <- MASS::stepAIC(lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log),
                        scale = 0,
                        direction = "both",
                        trace = 1,
                        keep = NULL,
                        steps = 1000,
                        k = 2)

summary(stepm2)
##============================================================================##
## Compare the two models
##============================================================================##
AIC(enter1, stepm1)
##============================================================================##
## A function for the relative importance
##============================================================================##
relweights <- function(fit,...){
R <- cor(fit$model)
nvar <- ncol(R)
rxx <- R[2:nvar, 2:nvar]
rxy <- R[2:nvar, 1]
svd <- eigen(rxx)
evec <- svd$vectors
ev <- svd$values
delta <- diag(sqrt(ev))
lambda <- evec %*% delta %*% t(evec)
lambdasq <- lambda ^ 2
beta <- solve(lambda) %*% rxy
rsquare <- colSums(beta ^ 2)
rawwgt <- lambdasq %*% beta ^ 2
import <- (rawwgt / rsquare) * 100
import <- as.data.frame(import)
row.names(import) <- names(fit$model[2:nvar])
names(import) <- "Weights"
import <- import[order(import),1, drop=FALSE]
dotchart(import$Weights, labels=row.names(import),
xlab="% of R-Square", pch=19,
main="Relative Importance of Predictor Variables",
sub=paste("Total R-Square=", round(rsquare, digits=3)),
...)
return(import)
}
##============================================================================##
##============================================================================##
## Check the relative importance
##============================================================================##
rel.1 <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, data = ryukyus_log)

relweights(rel.1, col="blue")
##============================================================================##
## The best subsets method
##============================================================================##
bes.tot <- leaps::regsubsets(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MaxD + MinD, 
                      data = ryukyus_log, 
                      nbest = 4)
##============================================================================##
##============================================================================##
## Visualise
##============================================================================##
plot(bes.tot, scale = "adjr2")
##============================================================================##
## Check the summary and the relative importance
##============================================================================##
tot <- lm(PD ~ Area + bio14 + bio8 + monthCountByTemp10 + MinD, data = ryukyus_log)  

summary(tot)


relweights(tot, col="blue")
sqrt(car::vif(tot)) > 2
##============================================================================##
## Compare the two models
##============================================================================##
AIC(enter2, stepm2, tot)
##============================================================================##
## Load the biotic cleaned data, as well as the spatial data for Italy                            
##============================================================================##
data_it_cleaned <- readRDS('Italian data cleaned.rds')

Italy <- read_sf('Italin provinces.shp') %>% 
  st_transform(4326) %>% 
  as_Spatial()
##============================================================================##
## First load the excel file with the correct species names
##============================================================================##
nms <- readxl::read_excel('italian_names.xlsx')
##============================================================================##


##============================================================================##
## Keep only one instance per taxon from our occurrence data
##============================================================================##
nms_flags <- data_it_cleaned %>% 
  distinct(scientificName) %>% 
  mutate(Taxon = nms$scientificName,
         action = nms$action) %>% 
  as_tibble() 
##============================================================================##

##============================================================================##
## First create the coords for the correct taxa
##============================================================================##
coords <- data_it_cleaned %>% 
  full_join(nms_flags) %>% 
  as_tibble() %>% 
  filter(!str_detect(action, 'DEL')) %>% 
  dplyr::select(decimalLongitude, 
                decimalLatitude,
                Taxon)
##============================================================================##


pt <- readRDS('Sparse matrix and relevant polygons for Italy.rds')
##============================================================================##
## Then create the sparse matrix
##============================================================================##
pt <- points2comm(dat = coords, 
                  mask = Greece, 
                  res = 0.04166667, 
                  lon = "decimalLongitude", 
                  lat = "decimalLatitude",
                  species = "Taxon")
head(pt[[1]][1:5, 1:5])
##============================================================================##
## Save it
##============================================================================##
saveRDS(pt, 'Sparse matrix and relevant polygons for Italy.rds')
italy_tree <- readRDS('Italian phylogenetic tree.rds')
##============================================================================##
## Create the phylogenetic tree
##============================================================================##
italy_tree <- data_it_cleaned %>% 
  full_join(nms_flags) %>% 
  as_tibble() %>% 
  filter(!str_detect(action, 'DEL')) %>% 
  dplyr::select(family, genus, Taxon) %>% 
  dplyr::rename(species = Taxon) %>% 
  distinct(species, .keep_all = T) %>%  
  get_tree(sp_list = .,
           tree = tree_plant_otl, # either 
           taxon = "plant", # or
           scenario = "S1",
           show_grafted = FALSE)


##============================================================================##
## Save it to disk
##============================================================================##
saveRDS(italy_tree, 'Italian phylogenetic tree.rds')
##============================================================================##
##============================================================================##
## Let's have a look on the p-values for each of the IVs that have entered
## our final model
##============================================================================##
summary(tot) 



##============================================================================##
## Inside the parenthesis we include the above p-values             
##============================================================================##
p <- c(0.0249, 0.0389, 0.2071, 0.2926, 0.0338) 
##============================================================================##


##============================================================================##
## Level of significance                             
##============================================================================##
alp <- 0.05 
##============================================================================##


##============================================================================##
## If the model's p-values are equal or smaller than the p-values we obtain 
## from the following command, then they are ok
##============================================================================##
p.adjust(p, method = "bonferroni", n = length(p)) 
##============================================================================##
## Create a new variable
##============================================================================##
sp_hot_cold_spots <- hot_cold_spots %>% 
  as_Spatial()
##============================================================================##


##============================================================================##
## Quick plot
##============================================================================##
plot(sp_hot_cold_spots,
     border = "grey",
     col = "lightgrey",
     main = "Diversity Hotspots and Coldspots")

plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$all_hot == 'HOT'), ],
     col = "red",
     add = TRUE, 
     border = NA)

plot(sp_hot_cold_spots[(sp_hot_cold_spots@data$all_hot == 'COLD'), ], 
     col = "blue", 
     add = TRUE, 
     border = NA)

legend("bottomleft",
       fill = c("red", "blue"),
       legend = c("hotspots", "coldspots"),
       bty = "n", 
       inset = .092)


##============================================================================##
## Second quick plot
##============================================================================##
hot_cold_spots %>% 
  ggplot() +
  geom_sf(aes(fill = cwe_sr_hot )) + 
  labs(title = 'Diversity hotspots in Italy',
       subtitle = 'Source: GBIF') + 
  scale_fill_manual(values = c('steelblue1','red', 'gray90'),
                    name = "Diversity hotspots",
                    labels = c("COLD", "HOT", "Neither"))
##===========================================================================##
## Do NOT run this code
##===========================================================================##
labs = knitr::all_labels(echo = TRUE)
##===========================================================================##

Βιβλιογραφία

Arponen, A. (2012) Prioritizing species for conservation planning. Biodiversity and Conservation, 21, 875–893.

Bottrill, M.C., Joseph, L.N., Carwardine, J., Bode, M., Cook, C., Game, E.T., Grantham, H., Kark, S., Linke, S., McDonald-Madden, E., & others (2008) Is conservation triage just smart decision making? Trends in ecology & evolution, 23, 649–654.

Buerki, S., Callmander, M.W., Bachman, S., Moat, J., Labat, J.-N., & Forest, F. (2015) Incorporating evolutionary history into conservation planning in biodiversity hotspots. Philosophical Transactions of the Royal Society B: Biological Sciences, 370, 20140014.

Cadotte, M.W. & Jonathan Davies, T. (2010) Rarest of the rare: Advances in combining evolutionary distinctiveness and scarcity to inform conservation at biogeographical scales. Diversity and Distributions, 16, 376–385.

Carta, A., Gargano, D., Rossi, G., Bacchetta, G., Fenu, G., Montagnani, C., Abeli, T., Peruzzi, L., & Orsenigo, S. (2019) Phylogenetically informed spatial planning as a tool to prioritise areas for threatened plant conservation within a mediterranean biodiversity hotspot. Science of The Total Environment, 665, 1046–1052.

Chape, S., Harrison, J., Spalding, M., & Lysenko, I. (2005) Measuring the extent and effectiveness of protected areas as an indicator for meeting global biodiversity targets. Philosophical Transactions of the Royal Society B: Biological Sciences, 360, 443–455.

Crisp, M.D., Laffan, S., Linder, H.P., & Monro, A. (2001) Endemism in the australian flora. Journal of Biogeography, 28, 183–198.

Daru, B.H., Roux, P.C. le, Gopalraj, J., Park, D.S., Holt, B.G., & Greve, M. (2019) Spatial overlaps between the global protected areas network and terrestrial hotspots of evolutionary diversity. Global Ecology and Biogeography, 28, 757–766.

Davies, T.J. & Yessoufou, K. (2013) Revisiting the impacts of non-random extinction on the tree-of-life. Biology letters, 9, 20130343.

Faith, D.P. (1992) Conservation evaluation and phylogenetic diversity. Biological conservation, 61, 1–10.

Farrier, D., Whelan, R., & Mooney, C. (2007) Threatened species listing as a trigger for conservation action. Environmental Science & Policy, 10, 219–229.

Forest, F., Grenyer, R., Rouget, M., Davies, T.J., Cowling, R.M., Faith, D.P., Balmford, A., Manning, J.C., Procheş, Ş., Bank, M. van der, & others (2007) Preserving the evolutionary potential of floras in biodiversity hotspots. Nature, 445, 757–760.

Good, R. (1845) The geography of the flowering plants. Longsmans,

Isaac, N.J., Turvey, S.T., Collen, B., Waterman, C., & Baillie, J.E. (2007) Mammals on the edge: Conservation priorities based on threat and phylogeny. PloS one, 2, e296.

Jackson, P.W. & Kennedy, K. (2009) The global strategy for plant conservation: A challenge and opportunity for the international community. Trends in plant science, 14, 578–580.

Jenkins, C.N. & Joppa, L. (2009) Expansion of the global terrestrial protected area system. Biological conservation, 142, 2166–2174.

Jetz, W., Thomas, G.H., Joy, J.B., Redding, D.W., Hartmann, K., & Mooers, A.O. (2014) Global distribution and conservation of evolutionary distinctness in birds. Current biology, 24, 919–930.

Kira, T. (1991) Forest ecosystems of east and southeast asia in a global perspective. Ecological Research, 6, 185–200.

Li, D., Trotta, L., Marx, H.E., Allen, J.M., Sun, M., Soltis, D.E., Soltis, P.S., Guralnick, R.P., & Baiser, B. (2019) For common community phylogenetic analyses, go ahead and use synthesis phylogenies. Ecology, 100, e02788.

Linder, H. (2001a) On areas of endemism, with an example from the african restionaceae. Systematic biology, 50, 892–912.

Linder, H. (2001b) Plant diversity and endemism in sub-saharan tropical africa. Journal of Biogeography, 28, 169–182.

Maekawa, F. (1974) Origin and characteristics of japan’s flora. The Flora and Vegetation of Japan., 33–86.

Margules, C.R. & Pressey, R.L. (2000) Systematic conservation planning. Nature, 405, 243–253.

Mishler, B.D., Knerr, N., González-Orozco, C.E., Thornhill, A.H., Laffan, S.W., & Miller, J.T. (2014) Phylogenetic measures of biodiversity and neo-and paleo-endemism in australian acacia. Nature Communications, 5, 1–10.

Nakamura, K., Suwa, R., Denda, T., & Yokota, M. (2009) Geohistorical and current environmental influences on floristic differentiation in the ryukyu archipelago, japan. Journal of Biogeography, 36, 919–928.

Ota, H. (1998) Geographic patterns of endemism and speciation in amphibians and reptiles of the ryukyu archipelago, japan, with special reference to their paleogeographical implications. Researches on Population Ecology, 40, 189–204.

Pollock, L.J., Thuiller, W., & Jetz, W. (2017) Large conservation gains possible for global biodiversity facets. Nature, 546, 141–144.

Reece, J.S. & Noss, R.F. (2014) Prioritizing species by conservation value and vulnerability: A new index applied to species threatened by sea-level rise and other risks in florida. Natural Areas Journal, 34, 31–45.

Rosauer, D., Laffan, S.W., Crisp, M.D., Donnellan, S.C., & Cook, L.G. (2009) Phylogenetic endemism: A new approach for identifying geographical concentrations of evolutionary history. Molecular ecology, 18, 4061–4072.

Smith, S.A. & Brown, J.W. (2018) Constructing a broadly inclusive seed plant phylogeny. American journal of botany, 105, 302–314.

Triantis, K.A., Guilhaumon, F., & Whittaker, R.J. (2012) The island species–area relationship: Biology and statistics. Journal of Biogeography, 39, 215–231.

Tucker, C.M., Aze, T., Cadotte, M.W., Cantalapiedra, J.L., Chisholm, C., Dı́az, S., Grenyer, R., Huang, D., Mazel, F., Pearse, W.D., & others (2019) Assessing the utility of conserving evolutionary history. Biological Reviews, 94, 1740–1760.

Tucker, C.M., Cadotte, M.W., Carvalho, S.B., Davies, T.J., Ferrier, S., Fritz, S.A., Grenyer, R., Helmus, M.R., Jin, L.S., Mooers, A.O., & others (2017) A guide to phylogenetic metrics for conservation, community ecology and macroecology. Biological Reviews, 92, 698–715.

Watson, J.E., Dudley, N., Segan, D.B., & Hockings, M. (2014) The performance and potential of protected areas. Nature, 515, 67–73.


  1. Βασική προϋπόθεση φυσικά είναι να τις έχετε εγκαταστήσει ήδη

  2. Πληκτρολογώντας vars στην κονσόλα σας

  3. Συγκρίνετε φερειπείν τα δεδομένα που έχουμε για το υψόμετρο με αυτά για οποιαδήποτε άλλη μεταβλητή

  4. Δεν λογαριθμούμε τη μεταβλητή PD. Γιατί;

  5. Στην προκειμένη περίπτωση αυτό έχει ήδη για εσάς - είναι μια αρκετά επίπονη διαδικασία

 

Delivered to you by Dr. Kostas Kougioumoutzis

kkougiou@aua.gr